%% 牛顿冷却定律温度求解案例:Newton cooling law temperature solution case
%% 数据初始化
k0=1;k1=0.99;k2=0.98;k3=0.9;k4=0.99;k5=0.01;% 传热系数
d=[0.6,0.6,3.6,0.6];% 厚度设置
d1=cumsum(d);
dx=[1e-2,1e-2,1e-1,1e-2];% 厚度微元设置
tw=75;% 外部温度75
tn=37;% 假人温度37
u1=tn*ones(1,length(0:dx(1):d1(1)));
u2=tn*ones(1,length(d1(1):dx(2):d1(2)));
u3=tn*ones(1,length(d1(2):dx(3):d1(3)));
u4=tn*ones(1,length(d1(3):dx(4):d1(4)));% 防护服温度初始化
tm=5400;% 总时间
u=zeros(1,tm);
for t=1:tm
    tic
    for i=1:length(0:dx(1):d1(1)) % 第一层防护服温度迭代
        if(i==1)
            deltau=-k0*(u1(i)-tw)-k1*(u1(i)-u1(i+1));
            u1(i)=u1(i)+deltau;
        elseif(i==length(0:dx(1):d1(1)))
            deltau=-k1*(u1(i)-u1(i-1))-k2*(u1(i)-u2(1));
            u1(i)=u1(i)+deltau;
        else
            deltau=-k1*(u1(i)-u1(i-1))-k1*(u1(i)-u1(i+1));
            u1(i)=u1(i)+deltau;
        end
    end
    for i=1:length(d1(1):dx(2):d1(2)) % 第二层防护服温度迭代
        if(i==1)
            deltau=-k1*(u2(i)-u1(end))-k2*(u2(i)-u2(i+1));
            u2(i)=u2(i)+deltau;
        elseif(i==length(d1(1):dx(2):d1(2)))
            deltau=-k2*(u2(i)-u2(i-1))-k3*(u2(i)-u3(1));
            u2(i)=u2(i)+deltau;
        else
            deltau=-k2*(u2(i)-u2(i-1))-k2*(u2(i)-u2(i+1));
            u2(i)=u2(i)+deltau;
        end
    end
    for i=1:length(d1(2):dx(3):d1(3)) % 第三层防护服温度迭代
        if(i==1)
            deltau=-k2*(u3(i)-u2(end))-k3*(u3(i)-u3(i+1));
            u3(i)=u3(i)+deltau;
        elseif(i==length(d1(2):dx(3):d1(3)))
            deltau=-k3*(u3(i)-u3(i-1))-k4*(u3(i)-u4(1));
            u3(i)=u3(i)+deltau;
        else
            deltau=-k3*(u3(i)-u3(i-1))-k3*(u3(i)-u3(i+1));
            u3(i)=u3(i)+deltau;
        end
    end
    for i=1:length(d1(3):dx(4):d1(4)) % 第四层防护服温度迭代
        if(i==1)
            deltau=-k3*(u4(i)-u3(end))-k4*(u4(i)-u4(i+1));
            u4(i)=u4(i)+deltau;
        elseif(i==length(d1(3):dx(4):d1(4)))
            deltau=-k4*(u4(i)-u4(i-1))-k5*(u4(i)-tn);
            u4(i)=u4(i)+deltau;
        else
            deltau=-k4*(u4(i)-u4(i-1))-k4*(u4(i)-u4(i+1));
            u4(i)=u4(i)+deltau;
        end
    end
    toc
    u(t)=u4(end);
end
